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Abstract 



Several integrate-to-threshold models with differing temporal in- 
tegration mechanisms have been proposed to describe the accumula- 
tion of sensory evidence to a prescribed level prior to motor response 
in perceptual decision-making tasks. An experiment and simulation 
studies have shown that the introduction of time-varying perturba- 
tions during integration may distinguish among some of these models. 
Here, we present computer simulations and mathematical proofs that 
provide more rigorous comparisons among one-dimensional stochastic 
differential equation models. Using two perturbation protocols and 
focusing on the resulting changes in the means and standard devia- 
tions of decision times, we show that, for high signal-to- noise ratios, 
drift-diffusion models with constant and time-varying drift rates can 
be distinguished from Ornstein-Uhlenbeck processes, but not neces- 
sarily from each other. The protocols can also distinguish stable from 
unstable Ornstein-Uhlenbeck processes, and we show that a nonlinear 
integrator can be distinguished from these linear models by changes in 
standard deviations. The protocols can be implemented in behavioral 
experiments. 



1 Introduction 



Rea ction t ime t asks have long been used to study human decision mak- 
ing (ILucd . Il986l ). One paradigm requires subjects to detect or discriminate 
sensory signals by making a voluntary motor response. Response time dis- 
tributions obtained in such perceptual-mo tor tasks allow inference of cogni- 
tive information processing fjPosnerl . Il978l ). With the advent of methods for 
recording neural activity in awake behaving animals, such tas ks have been 
adopted by neurophysiolo gists, especially in studies on monkeys (jSchalll . 12003 : 



Gold and Shadlenl . 120071 ). and neuronal firing rates in several brain areas 



have been shown to correlate with motor respons es. The lateral int r apari- 
etal area, frontal eye fields, and superior coUiculus (iGold and Shadlenl . 120071 ) 
all exhibit activities that ramp up over time toward a fixed level before a 
decision is signaled (e.g. b y a saccadic eye movernent in the direction of 
the r e corded response fi e ld) flHanes and Schalll . Il996l : iRoitman and Shadlenl . 
2002; Churchland et al. . 20081 ). The slopes of the ramps not only correlate 
with task difficulty (the harder the task, the lower the slope), but also with 
response time (higher slopes precede faster responses). These areas may 
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therefore provide neural substrates for integrating sensory information to- 
ward a decision criterion before a perceptual decision is made. 

Various integrate-to-threshold models have been p roposed to describe 
both response times and neurobiologi cal mechanisms fISmith and Ratclifl . 



20041). including drift- diffusion models (IRatclifj. Il978l: iMazurek et al 



2003; 



Smith and Ratclifj . l2004l : IPitterichl . l2006al : ISimen et al.l . l2006l:lRatcliff and McKoon 



2008 ) and attractor neural networks with mutual inhibition ( Brown and Holmesl. 



2001 lUsher and McClelland) . I2OOII : IWand . I2OO2I : IWong and Wand . l2006l : [Lo and Wand . 



20061 ). All share a common mechanism: accumulation or integration of sen- 



sory inputs toward a prescribed threshold, the first crossing of which deter- 
mines the decision and response time. These models may nonetheless be 
distinguished by the details of the integration process: drift-diffusion models 
typically accumulate evidence at a constant rate like a biased random walk, 
while attractor networks have unstable or stable steady- states that can, re- 
spect i vely, accelerate or decelerate the integration process (lUsher and McClellandl . 



2OOII : IWong and Wand . 120061 ). 



Often only behavioral data are collected for human subjects, and fits 
based on response times and choi ce accuracies a r e sornetimes unable t o distin- 
guish among competing models (IRatcliff et al.l . Il999l : iRatclifj . |2006| ). More- 
over, while cellular recordings in awake behaving animals offer direct insights 
into the integration process, fitting of both behavioral and neural data may 
still not suffice, especially when models incorporate multiple features and 
depend on multiple parameters. For examp l e, a dr i ft- diffusion model with a 



time-dependent 'urgency' signal (jPitterichl . l2006al : IChurchland et al.l . l2008l ) 



may be diff icult to separate from a recurrent netw ork model with strong self- 
excitation (IWand . I2OO2I : IWong and Wangl . |2006| ). Thus far, few principled 
attempts have been made to tease apart different integration mechanisms. 

Subthreshold electrical microstimulation of neural activities in beha ving 
animals may p r ovide more conclusive tests (jCohen and Newsomd . |2004| ). In 



Ditterch et al.l (120031 ). microstimulation of sensory neur ons was sho\yn to a f- 
fect the speed of the decision. More interestingly, in iHanks et al.l (120061 ). 
stimulation of cells with choice targets in the recorded neuronal response 
fields speeded decisions, but also reduced decision speeds when choice target 
directions were opposite to the recorded response fields, thus providing evi- 
dence of mutual inhibition. Pert urbations need no t be in vasive, and so can 
be used in human studies: e.g., iHuk and ShadlenI (120051 ) employed a brief 



motion pulse in the background of the primary visual stimulus; although it 
neither determined the choice nor influenced rewards, the pulse had a signif- 
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icant effect on response times. 



Earlier modeling efforts (IHuk and Sliadlenl.l2005l:IWong et al.l . 120071 : 1 Wong and Hukl . 

2008h have addressed the data of iHuk and ShadlenI (12005^ but they em- 



ployed many parameters and direct comparisons were not made between 
models. Here, we conduct a more rigorous study using simpler models, 
and seek more objective comparisons among them. We approximate four 
neural firing rate models as linear, scalar, stochastic differential equations 
(SDE), and we ask how their predicted response times are affected by short 
piecewise-constant perturbations with varying onset times and amplitudes. 
By comparing changes in means and standard deviations of response times, 
we demonstrate that the perturbations suffice to distinguish among the mod- 
els. Final ly, we show that a nonlin ear integrator model, qualitatively similar 
to that of I Wong and Wangl (120061 ). behaves much like one of the linear mod- 
els. 



2 Methods 



2.1 Reduction to one-dimensional linear integrate-to- 
threshold models 

Two-alternative forced-choice decision processes can, in essence, be modeled 
by two populations of excitatory neurons each endowed with self-excitatory 
connections and mutual inhibition via a shared inhibitory population. Each 
of the three populations can then be represented in a coarse-grained firing 
rate model by a single unit whose state describe s the population-av e raged 
activ ity rj{t) of the corresponding neuronal pool fl Wilson and Cowanl . 1 19721 . 
19731 ). Here we review the further reduction of a firing rate model, under 
suitable hy potheses, to a one-d i mensional dynamica l system. For additional 
details, see iBrown et all (120051 ): iBogacz et al.l (120061 ). 

Consider first the deterministic equations describing the firing rates of 
two excitatory (E) populations and a common inhibitory (J) population: 



te 



dri 

~dt 

dr2 

~dt 
dri 



-ri + FE(ri,r2,r/,/i) 
-rs + FE{r2,ri,ri,l2] 
-rj + Fi(ri,r2). 



(1) 
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Here te and tj are the synaptic time constants for the E and / units, 
and ri, r2 and rj are respectively the activity of neural units 1, 2 and /. The 
excitatory populations 1 and 2 are selective to stimuli 1 and 2, respectively. 
Fe and Fj are their input-output functions, and the overall input to each 
unit, with i = E or I, is li = Irecurrent,i + lstimuius,i- The decision time is 
the first passage time from stimulus onset to the first of ri or r2 reaching 
a prescribed decision threshold, which thereby signals choice 1 or 2. Since 
non-decision latencies (e.g. signal transduction and motor preparation) are 
usually assumed to be independent of stimulus strength I stimulus, we shall 
model the response time as the decision time plus a constant latency. Hence- 
forth we use the phrases 'decision time', 'reaction time' and 'response time' 
interchangeably. 

If the decision dyn amics passes near a sadd l e point and moves along 



Bogacz et al.l . l2006l ). the 



its unstable manifold (iBrown and Holmesl . 12001 
functions Fe, Fj may be linearized such that Eq. [1] simplifies to 

Te— = -ri + {ari- Pri + Ii + Iq^e), 
dt 



Te- 



dr2 

'dt 
drj 



-r2 + (ar2 - f3ri + I2 + /o,_ 



(2) 



Ti-^ = -ri + (7ri + + Jo,/). 



Here a,/?, 7 are the recurrent synaptic coupling strengths for self-excitation, 
inhibitory-to-excitatory, and excitatory-to-inhibitory connections, and Io,e 
(Iqj) is the constant background input to all the exc itatory (inhibitory) cells 
from outside the local circuit. Unlike IWangI (120021 ) . we exclude excitatory 
connections between ri and r2 and self-inhibitory connections, retaining only 
essential features. If tj <^te, we can further assume that the relatively fast 
dynamics of r/ equilibrates rapidly, such that ri ^ '~^{ri+r2) + Iq,i, and Eq. [2] 
becomes 



dv'i 

Te— = -n + (a - /37)ri - f3-fr2 + h + {Io,e - PIoj) 



Te- 



dr2 



-r2 + (a - (3-i)r2 - p-fn + h + (/o,£ - pIo,i) 



(3) 
(4) 



Defining a new variable X = ri—r2 and subtracting Eq. |4]from[3], we obtain 

dX 

— = kX + {h{t)-l2{t))/TE, 
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where k = {a — 1)/te contains the excitatory couphng strength and leak, 
and (Ji — I2)/te is proportional to the difference in inputs. The background 
inputs lo^E and Jq,/ and coupling strengths P and 7 cancel out. 

Generalizing the inputs Ji, I2 to be time- varying and including additive 
noise, the reduced dynamics is described by a one-dimensional SDE of the 
form 

dXt = b{Xt,t)dt + adWt, (5) 

where b{Xt,t) = kXt + (hit) — l2(t))/TE and a is the standard deviation 
of the noise, which is assumed to be a Wiener process with increments dWt 
drawn from a normal distribution with zero mean and unit variance. Eq. [5] 
provides a general description for noisy accumulator models, as illustrated 
schematically in Fig. [T]A.. 



2.2 Linear integrate-to-threshold models 

We are now in a position to introduce four linear instances of Eq. [51 When 
A; < (e.g. self-excitatio n a is relatively weak), we have a classical Ornstein- 
Uhlenbeck (OU) process (jUhlenbeck and Ornsteiru . ll930l : IWang and Uhlenbeck 
I945I ). If the decision threshold does not intervene, solutions eventually decel- 
erate and approach a stable steady-state X steady-state = {hit)—l2(t))/{\k\TE)- 
We denote this stable OU process SOU; in contrast, if /c > (stronger self- 
excitation), sample paths accelerate away from a fixed point in an unstable 
OU process, denoted UOU. SOU an d UOU provide reduced descriptions of 
the leaky competing accumulator of lUsher and McClelland! (120011 ) that re- 
spectively produce recency and primacy effects. 

In the special case of a = 1 (i.e. k = 0) an d constant stirn uli Ji, I2, Eq. [5] 
becomes a pure drift-diffusion (CD) equation fiRatchfj . Il978l ) 



dXt = bodt + adWt 



(6) 



in which the constant drift rate b = Bq = {Ii — I2)/te is proportional to the 
difference in stimuli. This CD model is simila r to, but differs in de tail from, 
the balance between leakage and inhibition in iBogacz et al. fl2006h .) 

In some variants of the CD model, drift rates can vary with time. This 
could be du e to dynamic stirn ul i, to some form of 'urgen cy' induced by the 
task des ign (Ditterich . 2006a| jbl: IChurchland et al" . 2008 ). to fluctuating at - 



tention (jSmith et al 



2004h or to other 'top-down' effects flLiu et al.l . |2008[ ) 
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Figure 1: One- dimensional integrate-to-threshold models and perturbation 
protocols. A: a general model for a two-alternative forced-choice reaction 
time task. Here the noisy state Xf, which quantifies the difference in firing 
rates between two competing excitatory cell populations, reaches the upper 
decision threshold in a sample trial, and choice is correct since threshold is 
in same direction as drift rate. Multiple sampling provides the response time 
distribution. B: the SPB model with unbiased input exhibits a subcritical 
pitchfork bifurcation at 5 = 0. Black curves show branches of stable (solid) 
and unstable (dashed) fixed points in the noiseless limit; branches for biased 
input shown in gray. C: a short pulse can be repeated over a range of onset 
times T after primary stimulus onset at t = to test effects on decision times. 
D: a pulse-antipulse perturbation whose leading pulse amplitude (red) can be 
adjusted to determine when decision times are unaffected by the combined 
pulse, yielding a zero-effect perturbation (ZEP). See text for further details. 



This time-dependent (TD) diffusion model is described by: 

dXt = b{t)dt + adWt. (7) 

We focus on the simple case b{t) = hot, since similar results obtain for other 
functions of time, but in Appendix B, we treat a modified TD model with 
time-varying perturbation amplitude. 
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2.3 A nonlinear integrate-to-threshold model 



The fi n al model to be consider ed is qualitatively similar to that of [Wang 
(I2OO2I ): IWong and Wang (120061 ) . It captures st ochastic dynamics in the 



neighborhood of a subcritical pitchfork bifurcatio n (iGuckenheimer and Holmes 
19831 : iBrown and Holmesl. I2OOII: IStrogata . I2OOII) that occurs as the stim ulus 
input level is varied (IRoxin and Ledberd . l2008l : IWong and Wand . |2006| ) : 



TxdXt = [eXt + Xf - Xf + b{t)]dt + adWf 



(8) 



Here h{t) = bo (constant) and e represent the biased and non-biased stimu- 
lus inputs respectively (Fig. [T]B). tx controls the overall temporal dynamics . 
Eg. p may be derived by norm al form theory (IGuckenheimer and Holmesl . 
1983 : Roxin and Ledbera . 20081 ): see Fig.[Tj3 for an illustration of the branches 
of stable and unstable fixed points in the noise- free limit. If |e| and \ b\ are 
suffici e ntly small, this rnqdel w i thout the Xf term approa ches that of [Wang 



(l2002h : IWong and Wand (120061 ) : lRoxin and Ledberd (120081 ). We include -X, 



to prevent unrealistic runaway activity, and, for simplicity, first assume a 
symmetrical system with 60 = 0, and then consider biased stimuli |6o| > 0. 
This model is denoted SPB. 



2.4 Perturbation protocols 

2.4.1 Single pulse perturbation with varying onset time 

We consider additive perturbations bi(t), under which Eq. [S] becomes dXt = 
[b{Xt,t) + bi(t)]dt + adWt, where bi(t) is applied from time T to T + AT 
and bi{t) = otherwise. We er nploy two piecewis e -cons tant forms, the first 



being the step function used by iHuk and ShadlenI (120051 ) 



0, t<T, 

biit) = { p, T <t<T + AT, (9) 
0, t > T + AT, 

in which the amplitude p can be positive or negative, assisting or opposing the 
unperturbed drift rate 60 due to the 'primary' stimulus. We fix the duration 
AT at 10% of the unperturbed mean first passage time tq/IO, apply bi at 
different onset times T and ask how the mean and standard deviation of the 
decision time change. 
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2.4.2 Zero-effect pulse-antipulse perturbation 



The second protocol uses a double pulse of the form 



hit) 



0, t < T, 

-Xp, T <t <T + AT/2, 

p, T + AT/2 <t<T- 

0, t > T + AT, 



AT, 



(10) 



where A is the relative height of the first pulse to the second. The second 
opposite-signed pulse attempts to reduce or cancel the expected change due 
to the first, and if the overall perturbation bi leaves the mean first passage 
time unchanged, we call it a zero- effect per turbation or Z EP. T his protocol, 
a variant of the paired pulse suggested in IWong et al.l (120071 ). is only ap- 
plied early in the integration process to avoid interference from the decision 
threshold. We shall seek critical values of A for which ZEPs occur for each 
model. Figs. [T]C and D illustrate both protocols. 



2.5 Simulations and parameter values 

Except for the SPB model, we first consider integration to a single decision 
threshold in the direction of the drift. This reduces the number of parameters, 
simplifies mathematical analyses, and helps isolate key effects. It applies to 
easy tasks in which drift rates are (relatively) high and errors rare. We 
subsequently relax this condition in our simulations. For the SPB model we 
first set Ji = I2, corresponding to difficult tasks and high error rates, and 
employ two thresholds. We then consider Ji 7^ I2 and a single threshold, 
representing easy tasks. 

Signal-to-noise ratios were chosen such that the variance in first passage 
(decision) times is significant, but not so great that an unreasonably large 
number of trials is needed to average out the noise, and we selected pertur- 
bation amplitudes, durations and onset times such that both protocols cause 
small but significant effects (e.g. maximum changes of ~ 10% in mean deci- 
sion times). In particular, d urations were an order of magnitude smaller than 



mean decision times Tq, as in lHuk and ShadlenI (120051 ). For the first protocol. 



onset times T were varied from stimulus onset at t = until there were no 
significant effects on mean decision time. In the second protocol we require 
T <ti tq. We set Xo = at t = to represent unbiased initial conditions, 
and the remaining parameters are chosen to ensure 'realistic' behavior (e.g. 
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the stable fixed point is above threshold for SOU; the unstable fixed point is 
below Xq = for UOU). Table [T] lists parameter values. 



Model 


b{Xt,t) 


Parameters 


c 


AT 


P 


CD 


bo + bi 


60 = 5, (T = 2.449 


20 


0.4 


5 


TD 


bot + bi 


60 = 4,(7 = 2.828 


20 


0.1 


4 


SOU 


kXt + bo 


k = -l,bo = 8, 


7 


0.4 


2 




+bi 


a = 1.414 








UOU 


kXt + 60 


k = 0.2, bo = 5, 


20 


1 


2 




+61 


a = 1.414 








SPB 


eX, + Xf 


e = -0.3, 0.05, feo = 


±0.75 


5 


0.005 




-XI + bo 


or 0.004, a = 0.01, 






or 0.0008 




+bi 


= 20 






{bo = 0.004) 



Table 1: Parameters used in simulations. Drift-diffusion model with constant 
drift rate (CD), and with time- varying drift rate (TD); SOU (UOU) stable 
(unstable) Ornstein-Uhlenbeck processes; SPB: nonlinear model. See Eqs. |5]- 
[8] and protocols [9] and [10] for details. Here ( is the distance from starting point 
X = to decision threshold, and for SPB e = —0.3 before stimulus appears 
at t = and e = 0.05 for t > 0. Perturbation durations AT and amplitudes 
p are values used in the first perturbation protocol. Parameters for the linear 
models with smaller signal-to-noise ratios are specified in Section 13.31 



We use a forward Euler-Maruyama scheme (IHighaml . I2OOII ) , integrating a 
sample path until it hits the prescribed decision threshold and recording the 
corresponding first passage time, at which the trial ends and the next begins. 
After collecting an appropriately large ensemble, we extract the first and 
second moments of the first passage time from the samples. The simulation 
is run once for each one of the different perturbations. For the models CD, 
TD, SOU and UOU, the step size is 10"'^ of a time unit and the sample size 
is = 10^, so that errors in the moments are of order 1/\/N ^ 10~^. Since 
changes in moments due to the perturbations are of order 10~^, these choices 
of time step and sample size reliably capture the effects of perturbation. The 
mean exit time for the unperturbed SPB model is ~ 80, so that a time step of 
0.01 and a sample size is 10^ suffices. We ran noisy simulations of this model 
with e = —0.3 for t < to represent pre-stimulus activity, and then switched 
to e = 0.05 for t > so that X = undergoes a pitchfork bifurcation and 
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becomes unstable, thus forcing a choice (see Fig. [TJ3. Magnitudes of e and 
a for SPB are much smaller than the corresponding feg and a for the linear 
models because the +Xf term accelerates solutions toward the thresholds). 

The parameters are not independent, and their number can be reduced by 
one in all models by dividing Eq. [5] by 60, and rescaling a o jh^, p p/bo, 
C C/bf), k k/bo and the dynamical variable Xt —>■ Xt/b^. The resulting 
CD and TD models are described by 2 parameters while the SOU and UOU 
models require 3 and the SPB model 4. Two additional parameters describe 
the amplitude p and duration AT of the perturbation. 

We focus on qualitative patterns of changes in response times and, to 
make unbiased comparisons, we adopt a dimensionless measure of the changes 
in means and standard deviations by normalizing the relative changes with 
respect to their unperturbed values. 



3 Results 

3.1 Single pulse perturbations of linear models 

Using a single pulse perturbation with varying onset time T, we estimated 
changes in mean decision time for the four linear integrate-to-threshold mod- 
els as described above, obtaining the results shown in Fig. [2l Imposed early 
during the integration process, a brief pulse in the same (opposite) direction 
as the drift rate significantly decreases (increases) mean decision times in all 
four models, as shown by the dashed (solid) curves. These effects fade as 
T increases due to a thresholding effect: for later perturbations, more tri- 
als have already crossed threshold, and so there are relatively fewer trials 
that are perturbed than unperturbed. Thus, when averaging over all trials, 
the infl uence of perturbat i ons is progressively reduc ed. Note that unlike the 



work of iHuk and Shadleru (120051 ) : IWong et all (120071 ). trials that have crossed 



threshold before and during perturbation are not excluded in the averaging 
process. This helps to reduce noise in the data, especially with late pertur- 
bation onset times. 

In Appendix A, we prove that this basic pattern - positive pulses advance 
mean decision times and negative pulses retard them - must hold provided 
that the perturbation occurs sufficiently early. The proof applies to a broad 
class of nonlinear systems including the SPB model of Eq. [HI provided that 
e > and thresholds lie inside the region in which drift magnitude increases 
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T/To 



Figure 2: Changes in mean decision times under single pulse perturbations 
with varying onset times superimposed on positive drift rates for linear in- 
tegrators CD, TD, SOU, and UOU (panels A-D). Vertical axes: normalized 
changes in mean first passage times (r — tq)/tq\ horizontal axes: normalized 
perturbation onset times T/tq. Changes due to negative pulses shown solid 
(— ); changes due to positive pulses dashed (+). CD and TD show near- 
constant changes for T/tq < 0.5 that are substantially greater for CD (note 
vertical scale on panel B), SOU exhibits a maximum at T/tq ~ 0.5 ("opti- 
mal" perturbation), and UOU shows monotonic decrease. Effects decrease 
in all CcLSGS cLS T approaches and passes Tq due to thresholding (see text). 



with X (i.e., well below the activation levels of the stable fixed points). 

Note that the drift-diffusion models (CD, TD) and the OU models (SOU, 
UOU) respond distinctly to early perturbations: the former exhibiting al- 
most constant changes in mean decision times (panels A and B), and the 
latter respectively showing increasing and decreasing effects (panels C and 
D). Moreover, due to thresholding, SOU is alone in having an onset time for 
which the perturbative effect is maximized. 

We next investigate how standard deviations are affected by perturbation 
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Figure 3: Changes in standard deviations of decision times under single pulse 
perturbations with varying onset times. All cases exhibit a maximum in 
normalized standard deviation (std(r)— std(ro))/std(ro), but only UOU has 
a (shallow) minimum. Effects of early pulses on all four integrators are 
distinct, as described in text. Format and conventions are as in Fig. O 



onset time: see Fig [31 Richer patterns arise than those of Fig. [2j Both 
CD and TD exhibit similar near-constant changes for early onsets, but the 
directions of changes for TD are opposite to those for CD, and the curves 
cross before reaching maxima. Standard deviations for SOU have a similar 
pattern to its mean, but UOU exhibits initial decreases followed by increases 
to a peak, unlike its mean. Optimal onset times exist for which the four 
standard deviations are maximally affected, and all four cases show distinct 
patterns. 

In the TD model considered above we assumed that the perturbation 
is not affected by the time-dependent gain b{t). This may not hold if the 
perturbation enters via the same sensory pathway as the stimuli. In Appendix 
B we consider a case in which perturbation and drift are affected in the same 
way. We show that this modified TD model differs qualitatively only for 
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early perturbation onset, for which it yields reduced effects. 



3.2 Zero-effect perturbations of linear models 

We now consider the zero-effect protocol of Eq. [TOl appealing to the fact 
that, in the noise-free case, the first passage time is a monotone function of 
the ratio A of the amplitudes of the opposing pulses. The critical A = A* that 
determines a ZEP is therefore given by the unique zero of 

r(A) = To, (11) 

where r(A) describes the functional dependence of passage time on A. Mean 
first passage times appear to remain monotone with respect to A when noise 
is included. Asymptotic methods can be u s ed to approximate passage time 
moments for small perturbations ( Lindner . 20041 ). but the lack of compact 



formulae for passage time distributions (PDFs) i n most cases, in cluding OU 



processes, and the need for double integrations (ILindnerl . |200J) makes this 
approach generally intractable. 

However, ZEPs may be analytically approximated for early onset times 
T with T' := T + AT ^ tq, such that the probability of threshold crossing 
prior to T' is negligible. In this sufficient condition is that the PDFs 

of the perturbed and unperturbed processes are identical at T', since for 
t > T', both processes are governed by the same drift and noise. They are 
therefore indistinguishable in the limit T, AT 0. Before deriving explicit 
expressions for A* under this assumption we present the results of numerical 
simulations with an intuitive explanation. 

Since the unperturbed SDE is independent of the current state Xf in 
the CD and TD cases, we expect that antisymmetric pulses with A* = 1, in 
which the pulses precisely cancel, will produce ZEPs. In contrast, responses 
to inputs decay with time for SOU, implying that the first pulse must be 
larger than the second for their net effect to cancel at perturbation offset 
t = T' (A* > 1). For UOU the reverse should hold (A* < 1). Fig. H confirms 
that this is the case. Table [2] lists the parameters used in the simulations; 
note that the early onset time conditions are only weakly satisfied (T' ranges 
from 25 -40% of tq). 

To derive explicit approximations for A*, we use a comparison method 
similar to that of Appendix A. Let Xt denote the unperturbed process and 
Zt the perturbed process. If Xq = Zq at t = 0, then Xt = Zt at the 
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Figure 4: Zero-effect perturbations of linear integrate-to-threshold models 
CD, TD, SOU, UOU (panels A-D). Vertical axes: normalized change in mean 
first passage times (r — ro)/ro; horizontal axes: relative amplitude A of first 
pulse compared to second. Zeroes of solid lines (— , positive pulse followed by 
negative pulse, p < 0) and dashed lines (+, negative pulse followed by positive 
pulse, p>0) identify ZEPs at A* ^ 0.9987, 1.0039 (CD); 0.9952, 1.0004 (TD); 
1.2271,1.2230 (SOU) and 0.9066,0.9034 (UOU). 



onset time T. When the perturbation ends at t = T' the solution of the 
unperturbed SDE of Eq. [S]with constant drift rate bo is 

pT+AT 

Xt' = {Xt + ho/k)e''^^ - h^/k + a e'^^^'-'UWs, (12) 

Jt 

while the perturbed system satisfies 

„T+AT rT+AT 

Zt, = {XT+bo/k)e''^^-bo/k+ / e''^'^^~%i{s)ds+a / e^^'^'-'^dW,. 

Jt Jt 

(13) 

These expressions well approximate the true activity levels for models CD, 
SOU and UOU with absorbing thresholds only when both processes have 
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Model 


To 


std(ro) 


T 


AT 


P 


A* (thy.) 


A* (sim.) 


CD 


4.005 


0.960 


0.5 


0.5 


5 


1 


1.0013 


TD 


3.145 


0.420 


0.5 


0.5 


5 


1 


0.9978 


SOU 


1.831 


0.374 


0.1 


0.4 


2 


1.2214 


1.2251 


UOU 


2.954 


0.142 


0.2 


1.0 


2 


0.9048 


0.9050 



Table 2: Parameters used in ZEP simulations, tq and std(ro) are the first 
and second moments of the passage time without perturbation, T and AT 
are the onset time and total duration of the perturbation; p is the amplitude 
of the second pulse and A* (thy.) and A* (sim.) are the ZEP pulse ratios 
predicted by the theory and obtained by averaging p < and p > results 
of Fig. H 



low probability of hitting the threshold during the interval [0,T'], but this 
condition holds for T'/tq ^ 1, since the integrated drift term is almost zero 
at early times. Eqs. [T2] and [T3] differ in the term e^^"^ e~^^bi{t)dt, which 
enters Zt' due to the action of the perturbation bi during the interval [T, T']. 
The ZEP condition for arbitrary perturbations is therefore 

T+AT 

e-%i{t)dt = 0, (14) 
which for the piecewise constant pulses of Eq. ITU] implies that: 

Hence for A; > (UOU), A* < 1, while for A; < (SOU), A* > 1. In the 
special case of CD, k = and A* = 1. 

For TD the drift term is time- dependent but the linear term kXt is absent 
and so the corresponding solutions of the SDEs are 

Xt' = / bo{s)ds + a dWs 
Jt Jt 

and 

pT+AT pT+AT pT+AT 

Zt' = / bo{s)ds+ I bi{s)ds + a j dW^. 

JT JT JT 
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In this ZEP must satisfy 

cT+AT 



hi{t)dt = 0, 



(16) 



T 



which yields A* = 1 for the pulses of Eq. (TDl as for CD. This result h o lds fo r 
all time-depende n t drift rates, including t hose used in ISmith et al.l (l2004j ) ; 
DitterichI fl2006bh : Iciuirchland et all fl2008h . 



Summarizing, ZEPs occur for the CD and TD models when the opposing 
pulse amplitudes are equal (A* = 1), but for SOU and UOU A* = e"'"^^/^ is 
larger and smaller than one, respectively. The ratios A* predicted by Eq. [TS] 
for the parameters of Table [2] are given in the final column of the Table. In 
all four cases they agree well with the zero crossings of the mean first passage 
times in Fig. |U in spite of the fact that T'/tq ^ 0.25 — 0.4 is not very small. 



3.3 Signal-to-noise ratio can influence perturbation ef- 
fects 

To what extent do the results of Section 13.11 hold when signal-to-noise ratios 
are reduced and a second threshold is added to track errors? To investigate 
this question we select new parameter values according to the criteria of 
Section 12. 5[ 

For the CD model, we increase a from 2.449 to 2.828 and set thresholds 
at C = ±5 instead of ±20, yielding an error rate of 5%. This introduces more 
complex behavior in which correct and error trials respond to the pertur- 
bations in differ ent manners, as shown i n Fig . [5]A. The combine d averages 



as computed in iHuk and ShadlenI (120051 ) and IWong et al.l (120071 ). see black 



solid and dashed curves) preserve some features of previous results: for early 
perturbations the direction of changes in mean RT agrees with Fig. [2JA., but 
the approximately-constant change in mean decision time early in the trial 
is replaced by a monotonic decline. Increased noise yields earlier threshold 
crossings, advancing the thresholding effect and masking the signatures of 
Fig.[2JA.. Moreover, the change in mean error RTs can reverse sign for late per- 
turbation onsets. We checked this seemingly counterintuitive phenomenon 
numerically by fixing the random generating seed for noise in a sample trial, 
finding that perturbation of long RT trials can change an impending error to 
a correct choice (data not shown). 

Increasing the noise level to cr = 7.071 with thresholds ( = ±20 produces 
an error rate of 10% in the TD model and changes in mean RTs for correct 
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T/Xo 



Figure 5: Changes in mean decision times for linear models with lower signal- 
to-noise ratios. Perturbation protocol is as in Fig. [2l Green (red) curves: 
perturbation p in the same (opposite) direction as drift rate Bq; solid: correct 
choices; dash-dotted: error choices. Solid (dashed) black curves are averages 
of correct and error choices with perturbation in the same (opposite) direction 
as drift rate. Grey (black) curves in panel A are for error rate of 0.02% (5%). 



and error choices similar to those of CD. However, averaged over correct and 
error choices, overall changes in mean RT exhibit a pattern similar to that 
of Fig. [213 (see black curves in Fig. [5|3)- In this respect TD maintains its 
signature more robustly with increasing error rate than CD, although the 
effects are progressively masked. 

For the SOU model, we increase a to 6.325, yielding an error rate of 
17%. The optimal perturbation onset time for changes in mean RT averaged 
over correct and error trials is masked (compare Fig. [2p with black curves 
in Fig. inP). This error rate is relatively high, but we have checked that for 
error rates as small as 3% SOU model features can also be masked, indicating 
that, like the CD model, SOU is sensitive to noise. 
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For the UOU model, more parameters must be adjusted to obtain rea- 
sonable RT changes. Here we set a = 2, ( = ±10, k = 0.02 and Bq = 0.5, 
yielding a high error rate of 23.5%, but similar trends occur for error rates 
from 0.1% to 35%. Fig. [5p shows that these are comparable to the previ- 
ous UOU results with high signal-to-noise ratio. Overall, we are unable to 
distinguish between noise-masking effects and intrinsic features of the UOU 
integrator. 

Using the same parameters, we find that changes in standard deviations 
for all the linear models show optimal perturbation onset times (Fig. [6]). The 
CD, SOU and UOU models are mutually indistinguishable (cf. Figs. [6]A., C 
and D), but TD exhibits a crossing effect (Fig. EP) similar to Fig. [313. 




T/Xn 



Figure 6: Changes in standard deviations of decision times for the linear 
models with lower signal-to-noise ratios. Format is as in Figs. [3] and [51 
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^ Change in mean g Change in standard 
o.4n decision time 0.2-1 deviation of decision time 



0.2 H 


-0.2 H 
-0.4 



I I I I I 

0.2 0.4 0.6 0.8 1 



0- 
-0.2- 
-0.4- 

-0.6- 



■ - +or + - 
+ + or 



I I I I I 
0.2 0.4 0.6 0.8 1 




0.4 0.6 0.8 1 1.2 1.4 1.6 0.4 0.6 0.8 1 1.2 1.4 1.6 

Figure 7: Perturbations of means (left column) and standard deviations 
(right column) of first passage times for the SPB model. Panels A, B: effects 
of single pulse perturbations with varying onset times; vertical axes as in 
Figs. M and [3l — h denotes effect of negative pulse at positive threshold, H — 
denotes effect of positive pulse at negative threshold, etc. Due to symme- 
try, — h and H — have same effect (solid curves) as do ++ and (dashed 

curves). Panels C, D: responses to zero-effect perturbations; vertical axis 
in C as in Fig. HI vertical axis in D shows normalized change in standard 
deviation as in Fig. [3l 



3.4 Perturbations of the nonlinear model 

We next apply both perturbation protocols to the nonlinear SPB model. The 
mean and standard deviation of decision times for the unperturbed system 
are 80 and 22 time units, respectively. The parameters used for the first 
protocol are given in Table [H while for the ZEP we set p = 0.005, T = 30 
and AT = 5, consistent with application early in the accumulation process. 
Figs. [7]\ and B show that the means and standard deviations of decision 
times behave like those for UOU (cf. panels D of Figs. [2] and [3]), decreasing 
monotonically and approaching zero as onset times increase. ZEPs occur in 
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the pulse ratio range 0.8 < A* < 1 for both mean and standard deviation 
(Figs. [7p and D). Assuming a hnear change in mean in Fig. [7p, we find that 
A* ^ 0.87. 




Figure 8: Changes in mean (A) and standard deviation (B) of decision times 
of the SPB model with high signal-to-noise ratio. All choices are correct. 
Perturbation protocols and plot formats are as in Figs. [7JA. and B. 



To compare with the behavior of the linear systems in the single thresh- 
old case, we now modify the SPB model so that it makes as few errors as 
they do. This is done by changing bo from to 0.004, which breaks the 
symmetry of the subcritical pitchfork bifurc ation (Fig. [JB) so th a t dec isions 
favor one choice, representing easy tasks (cf. iRoxin and LedbergI (120081 )): all 
other parameters remain the same. Fig. [8] shows that the basic behavior of 
Figs. [7p and D is preserved: effects on both mean and standard deviation of 
RTs decrease monotonically as perturbation onset time increases. While the 
effects on mean passage times of SPB are much like those for UOU, the fact 
that standard deviations of all four linear models exhibit maxima (Fig. [3]) 
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potentially allows us to distinguish SPB from them. The effects of ZEPs 
in this biased regime is similar to that shown in Figs. [7p and D (data not 
shown) . 



4 Discussion 



In this work, we have investigated how several simple stochastic integrate- 
to-threshold decision-making models respond to short perturbations. We 
focus on two-alternative forced-choice reaction time tasks for which general 
firing rate models with three coupled neural populations can be reduced 
to one-dimensional integrate-to-threshold systems. These reduced systems, 
which may be pure drift- diffusion or Ornstein-Uhlenbeck processes, can be 
compared analytically and in simulations that require few parameters. 

We examined two perturbation protocols, the first being a brief pulse with 
variable onset time, and the second a pulse-antipulse pair whose amplitude 
ratio can be adjusted to produce minimal effects (the zero-effect perturba- 
tion, ZEP) . The simulations of §3] (Figs. [2] and [3]) show that the changes 
in means and standard deviations of decision times allow both perturbation 
protocols to distinguish between drift-diffusion models and stable and unsta- 
ble Ornstein-Uhlenbeck processes, provided that perturbations are delivered 
sufficiently early in the integration process and the signal-to-noise ratio is not 
too small. Changes in standard deviations can also assist in distinguishing 
drift-diffusion models with constant drift from those with time-dependent 
drift rates (cf. Figs. [3]A. and B), alt hough diff e rences are small and may not 
be detectable experimentally. (As iDitterichI (j2006aj) notes, entire reaction 
time distributions can further assist in this.) For early perturbations the 
ZEP conditions can be approximated analytically, making predictions that 
are confirmed by computer simulations (Fig. Hj). While general analytical 
expressions for systems with finite noise appear elusive, in Appendix A we 
prove inequalities for more general nonlinear systems that partially explain 
our simulation results. 

Although we focus on linear models, our methods extend to other sys- 
tems, and we include results for a nonlinear model that capture the dynamics 
near a s ubcritical pit c hfork bifurcation which is typical of reduced popul ation 
models (jWangj . I2OO2I : IWong and Wand . I2OO6I : iRoxin and Ledbergj . boosi ). As 
the results of Appendix A suggest, it behaves in a manner similar to the un- 
stable OU model under both perturbation protocols, although perturbative 
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effects on standard deviations may allow distinctions to be made (Figs. [31H]). 
We have also confirmed that the TD model's behavior remains similar with 
other increasing drift rates (e.g. b{t) ~ bot"^; simulation r esults not sh o wn). 

These results reinforce the simulations and claims in IWong et al.l (120071 ) 
that time-varying perturbations can in principle reveal integration mech- 
anisms. However, some parameter ranges frustrate clear-cut comparisons 
among models. The first occurs for low signal-to-noise ratios and high er- 
ror rates (cf. Figs [MS]), and so can be mitigated by using perturbations 
with more easily discriminated stimuli. Secondly, if there is a 'dead time' 
between stimulus onset and the start of evidence integration, the CD and 
TD models can display similar signatures to the SOU model, even with high 
signal-to-noise ratios. 

Since perturbations can be delivered through the senses as well as by 
direct electrical stimulation, the method can be noninvasive and is there- 
fore appropriate for human subjects. Furthermore, unlike fitting neuronal 
firing rates, both p erturbation protocols can distinguish linear and nonlinear 



UOU - type models flWangj. 120021 : IWong and Wangl . 120061 : IRoxin and Ledbergl . 



20081 : IWong and Hukl . |2008| ) from dri ft-diffusion models with an ' urgent 



time-dependent increase in drift rate (jSmith et al.l . l2004l : iDitterichl . l2006a 



Churchland et al.l . l2008l ). However, care is required in designing such exper- 



iments. Although a high signal-to-noise ratio can in principle help identify 
integrators, the reaction times can become too short, masking the signature 
behaviors of the integrators with the 'thresholding effect'. A possible solu- 
tion could be to prolong the response deadline in relatively simple reaction 
time tasks and analyze only data from long RT trials. 

Perturbation inputs to decision-making circuits may originate in a single 
sensory pathway or in other br ain areas, as i n mu l tisensory and 'to p-down' 
inputs (e.g. due to attention fISmith et al.l . 12004 : iLiu et al.l . |2008|)) which 



could modulate a decision throughout a sensorimotor pathway. The general 
notion of perturbation developed in the reduced model context extends to 
such modulatory inputs. The methods may also g eneralize to other tasks 



that recruit neural i ntegra tors (I Goldman et al. 
dShea-Brown et al.l . [2006h . 



20081 ). such as interval timing 
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Appendix A: Rigorous estimates on first pas- 
sage times 

In this appendix we analyze short single pulse perturbations, First in the 
noiseless limit, and then with additive noise. The inequalities proved here 
suggest explanations for the monotonicity of mean first passage times of SOU 
and UOU processes with respect to early onset times (cf. Figs. [2]C and D); 
moreover, they apply to more general, nonlinear one-dimensional systems 
such as SPB with thresholds \(\ sufficiently close to Xt = 0. 

The noiseless limit 

As noise amplitude tends to zero, the SDK Eq. [5] becomes an ordinary dif- 
ferential equation (ODE): 

rl X 

^:=A- = /W. (17) 

Consider two compact perturbations p and p with the same amplitude profile 
g{t) and duration AT, applied at different onset times Ti and Ti: 

p{t) = l[nMt)9{t - Ti), Pit) = - Ti) (18) 

(here 1[Ti,T2] denotes the indicator function that takes the value 1 for Ti < 
t < T2 and zero otherwise). We choose < Ti < T2 < Ti < T2 and 
T2 — Ti = T2 — Ti = AT such that p and p do not overlap and that neither 
perturbed solution reaches the threshold X = ( before t = T2. We shall 
compare solutions Yf and Yf of the corresponding perturbed ODEs, 

Y = f{Y)+p{t) and Y = f{Y)+p{t), (19) 

started at the same initial condition Yq = Yq = 0, to show how the signs of 
/' and g determine their threshold passage times r and f . We allow smooth 
nonlinear functions /, but assume that f{X) > for all < X < C and that 
g does not change sign. The fact that solutions of scalar ODEs cannot cross 
each other is a key tool. 

Lemma 1: Let Xt{x) denote the solution of Eq. [7^ with initial value x. 
Then Xt{xi) < Xt{x2) for any t > if and only if xi < X2- 
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Proof: This follows from the uniqueness of solutions of ODEs and the or- 
dering of points in a one-dimensional phase space. □ 



Since neither perturbation acts after t = T2, and neither solution has 
reached threshold at T2, Lemma 1 implies that r < f if and only if Y{T2) > 
Y[T2). We now prove the main result summarized in Table [3l which is a 
corollary of Proposition 2. 





/' > L > 


f'<-L<0 


g>0 


T < T 


T > T 


g<0 


T > T 


T < T 



Table 3: Monotonicity of first passage dependence on onset time for four sign 
combinations of /' and g, in the noise-free limit. L is some constant. 

Proposition 2: Assume that the function g{t) in Eqs. [721 is strictly positive 
for t G (0, AT). Then for some constant L (a) If f > L > 0, then Y(t) > 
Y{t) for all t > T2; (b) If f < -L < 0, then Y{t) < Y{t) for all t > Ta- 

Proof of part (a): We need only establish the inequalities at t = T2, since 
no perturbations occur after this time. We compare solutions of Eqs. [19] 
using the fact that /(X) > and f'{X) > L > for X > implies that 
f{Y) - f{Y) > LiY - Y). Letting u{t) = Y{t) - Y{t) we have 

u = f{Y) - f{Y) + pit) - pit) >Lu+ Pit) - Pit). (20) 

Since the solutions coincide until the first perturbation occurs at t = Ti, 
m(Ti) = 0, and integrating inequality [20] over the interval [Ti,T2], on which 
pit) = git — Ti) and pit) = 0, yields: 

^(^2) > / ' e^^^'-'^git - Ti)dt > 0. (21) 

Both pit) and pit) = on [T2,Ti], and integration of [20] and use of [2T] gives 

m(Ti) > n(T2)e^(^^-^^) > e^^^^'^^) e^(^^-*)^?(t - T,)dt > 0. (22) 
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Finally, using the fact that p{t) = on [Ti,!^], we integrate [2U] again to 
obtain 



gL(f2-T2) _ 



rT2 



> 



AT 



g{t - fi)dt 



gL(AT-t)^(^)^^ > 0, 



(23) 



where we also use inequality [22] and the fact that the additive perturbations 
of Eq. [IHlare identical over the time intervals [Ti,T2] and [Ti,T2] to rewrite 
the integral term. Hence Y{T2) > Y{T2) and r < f. 

Proof of part (b): We again compare solutions of Eqs. [T2]but the direction 
of the ineq ualities is now reversed, so that we are in essence using Gronwall's 
inequality (IGuckenheimer and Holmesl . Il983l ). In this case Eq. [201 becomes 



ii < —Lu + pit) — pit), 



(24) 



inequality [22] becomes 

u{f,) < «(T2)e-^(^^-^^) < e-^(^^-^2) 
and the final inequality reads 



rT2 



-L{T2~t) 



git - Ti)dt < 0, (25) 



u 



(T2) = F(T2)-F(r2)< 



1 - e 



-L{f2-T2) 



e-^(^^-*)(7(t)t/t < 0. (26) 



This implies that r > f, as claimed. Proofs are similar for < 0. 



□ 



Extension to noisy systems 

The conclusions of Proposition 2 and Table [3] may be extended to apply to 
expected passage times for SDEs with additive noise. Specifically, replacing 
Eqns. [Hby 

dYt = [fiVt) + pit)]dt + adWt, (27a) 
dYt = [f{Yt) + p{t)]dt + adWt, (27b) 

we have the following result: 
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Proposition 3: Assume that the function g{t) in Eqs. [72| is strictly positive 
for t G (0, AT). Then for some constant L (a) If f > L > 0, then E[r] < 
E[f]; (b) If f < —L < 0, then E[r] > E[f], where r and f denote the first 
passage times for Eqns. \27a\ and 27b respectively. 

Proof: We first compare solutions Yt and Yt obtained from Eqns. I27al and 
12 7bl under tlie same sample noise path with increments dWt = dWt. Since the 
stochastic integral dWt is almost surely continuous ( Feller . 1957 ). and the 
integrated noise terms cancel, the comparison of solutions using u{t) = Yf — Yt 
proceeds as in the deterministic case above, and we conclude that t < f 
ioT g > and /' > L > 0, and r > f for 51 > and f < -L < 0. 
Finally, averaging each process over an ensemble of sample paths, we may 
conclude that the mean first passage times of the two processes satisfy similar 
inequalities to those in the top row of Table [31 E[r] < E[f] and E[r] > E[f] 
respectively. A similar argument applies to the case g < 0. □ 



Appendix B: A TD model with time- varying 
perturbation amplitude 

In the body of the paper, perturbations of the TD model were assumed 
unaffected by the time-dependent gain. Here we extend the simulations to a 
simple case in which perturbation and stimulus (drift) share the same overall 
linear increase in amplitude. Specifically, we assume the dynamics to be 
governed by 

dXt = (bo + p{t))tdt + adWf (28) 

Otherwise, we change only the perturbation amplitude, taking p = ±2 in- 
stead of 4. Fig. [9JA. shows that the changes in mean and standard deviations 
of exit time for this model remain similar to the previous TD model, but with 
smaller perturbative effects near the beginning of the trial (quite similar to a 
SOU model) due to the monotonic time-dependence of perturbation ampli- 
tude. However, although the previous pattern of mean exit times no longer 
applies, the behavior remains qualitatively dissimilar to the CD, UOU and 
SPB models. 
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Figure 9: Normalized changes in mean (A) and standard deviation (B) of 
first passage times for single pulse perturbations with varying onset times 
applied to a variant of the TD model. Vertical axes and curve conventions 
are as in Figs. [2] and O 
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